function [segI, loc] = detectBall(I)
    % load('mu.mat');
    load('Samples.mat');
    mu = mean(Samples); % mu
    load('sig.mat');
    
    Id=double(I); % array in size of (row, col, 3)
    
    row=size(Id,1); 
    col=size(Id,2);

    % x_i - mu
    for i=1:3
        muval = mu(i);
        Id(:,:,i) = Id(:,:,i) - muval;
    end
    
    % reshape the image to a matrix in size of (row*col, 3)
    Id=reshape(Id,row*col,3);
     
    % calc possibility using gaussian distribution
    % be careful of using * and .* in matrix multiply
    Id = exp(-0.5* sum(Id*inv(sig).*Id, 2)) ./ (2*pi)^1.5 ./ det(sig)^0.5;
    
    % reshape back, now each pixels is with the value of the possibility
    Id=reshape(Id,row,col);
    
    % set threshold
    thr=8e-07;
    
    % binary image about if each pixel 'is ball'
    Id=Id>thr;
    
    % find the biggest ball area
    segI = false(size(Id));
    
    CC = bwconncomp(Id);
    numPixels = cellfun(@numel,CC.PixelIdxList);
    [biggest,idx] = max(numPixels);
    segI(CC.PixelIdxList{idx}) = true;
    %figure, imshow(segI); hold on;
    
    S = regionprops(CC,'Centroid');
    loc = S(idx).Centroid;
    %plot(loc(1), loc(2),'r+');